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Abstract: The determination of parameters of controllers is an 
important problem in automatic control systems. In this paper, 
the Levenberg Marquardt (LM) Algorithm is used to effectively 
solve this problem with reasonable computational effort. The 
Levenberg Marquardt (LM) Algorithm for optimization of 
three term (PID) controller parameters with dynamic model of 
pH neutralization process is presented. The main goal is to 
show the merits of Levenberg Marquardt algorithm 
optimization and to determine its suitability in the area of 
control systems. Lastly, the application of this approach to the 
calculation of the parameters of PID controller shows that the 
Levenberg Marquardt (LM) algorithm has a better dynamic 
performance of pH neutralization process. 
Keywords: Parameters Optimization, PID controller, pH 
neutralization Process, Levenberg Marquardt (LM) 
algorithm 

1. Introduction 

In mathematics and computing, the Levenberg-Marquardt 
Algorithm (LMA), also known as the damped least-squares 
(DLS) method, is used to solve non-linear least squares 
problems. These minimization problems arise especially in least 
squares curve fitting. The LMA interpolates between the Gauss- 
Newton Algorithm (GNA) and the method of gradient descent. 

The LMA is more robust than the GNA, which means that 
in many cases it finds a solution even if it starts very far off the 
final minimum. For well-behaved functions and reasonable 
starting parameters, the LMA tends to be a bit slower than the 
GNA. LMA can also be viewed as Gauss-Newton using a trust 
region approach. The LMA is a very popular curve-fitting 
algorithm used in many software applications for solving generic 
curve-fitting problems. However, as for many fitting algorithms, 
the LMA finds only a local minimum, which is not necessarily 
the global minimum [ 1 ] . The primary application of the 
Levenberg-Marquardt algorithm is in the least squares curve 
fitting problem: given a set of m empirical datum pairs of 
independent and dependent variables, (xi,yi) optimize the 
parameters /? of the model curve /(*,/?) so that the sum of the 
squares of the deviations becomes minimal. 

m 

i=i 



The Levenberg-Marquardt (LM) algorithm is an iterative 
technique that locates the minimum of a multivariate function 
that is expressed as the sum of squares of non-linear real-valued 
functions [4, 6]. 

In the other hand, a common problem in control system 
design is establishing the appropriate value of controller gains. In 
general a low value of gain produces a slow system response, 
while high gain values can cause an excessively-oscillatory 
response with the possibility of instability. Somewhere between 
these extremes is a value of gain that produces the best system 
response. The essential function of a feedback control system is 
to reduce the error, e(t) between any variable and its demanded 
value to zero as quickly as possible. Therefore, any criterion 
used to measure the quality of system response must take into 
account the variation of the error over the whole range of time. 
Four basic criteria are in common use: 

Integral of absolute error (IAE) = j^\e(t)\dt 
Integral of squared error (ISE) = / Q c {e(t)} 2 dt 
Integral of time multiplied by absolute error 

(ITAE) = j™t\e(t)\dt 

Integral of time multiplied by squared error (ISE) 

= ^t{e(t)fdt 

For any of the possible criteria, the best response 
corresponds to the minimum value of the chosen criterion. Note 
that in all cases it is either the absolute error or the squared error 
which is involved, straightforward integration of the error would 
produce zero result, even if the system response was a constant 
amplitude oscillation. IAE is often used where digital simulation 
of a system is being employed, but it is inapplicable for 
analytical work, because the absolute value of an error function 
is not generally analytic in form. This problem is overcome by 
the ISE criterion. The ITAE and ITSE have an additional time 
multiplier of the error function, which emphasizes long-duration 
errors, and therefore these criteria are most often applied in 
systems requiring a fast settling time [3, 5]. 

In this paper, the shape of the complete closed loop 
response, from time t — 0, until steady state has been reached, 
could be used for the formulation of a dynamic performance 
criterion. The simple criteria of this category are based on the 
entire response of the process and the integral of the Square 
Error (ISE) criterion used here, where 

p CO pCG 

ISE= e 2 (t)dt= {Y sp (t)- Y(t)} 2 dt (1) 

■>0 ^0 
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Where e(t) = Y sp (t)- Y(t) is the deviation (error) of the 
response from the desired set point. 

The ideal continuous time domain PID controller for a 
SISO process is expressed in the Laplace domain as follows: 



K, 



G(s) = K p + — + K d S 



(2) 



With K p = proportional gain, K, = integral time constant and K d 
= derivative time constant. 

2. Mathematical model of pH neutralization process 

Consider a pH neutralization process as shown in Fig. 1 . The 
flow rates of acid, buffer, base and effluent streams are denoted 
by (ft, 42»?3» and 174, respectively. Output of the process is the 
pH value of the effluent stream, and the flow rate of base 
stream, q 3 is the control input. A dynamic model is derived using 
the conservation laws and reactions equilibrium. The modeling 
assumptions include perfect mixing, constant volume of the 
neutralization tank (V), and complete solubility of the ions 
involved. The chemical reactions in the system are as follows 
[8]: 



q 2 (buffer) :NaHC0 3 



q,(acid) :HN0 3 +H 2 CO 



q 3 (base) :NaOH+NaHCO; 





q /effluent) 



Fig. 1 pH neutralization process 
H 2 C0 3 ^ H + + HC0 3 
HC0 3 <-» H+ + C0\~ 
HN0 3 ->//+ + JVO3- 
NaHC0 3 -» Na+ + HC0 3 
NaOH -> Na + + 0H~ 

The equilibrium constants for these reactions are: 
[HC0 3 ][H + ] 

KCLi 



Ka 2 = 



[H 2 C0 3 ] 



\HC0o 



Kw = [H + ] + [0H~] + [C0|-] 

The chemical equilibrium equations are modeled using the 
reaction invariant concept. For this system, concentrations of 
reaction invariants are defined as: 



*i = [NO3] 
x 2 = [Na + ] 

x 3 



[H 2 C0 3 ] + [HC0 3 ] + [COt] 



(3) 



Denoting, y — pH the ions neutrality balance in the tank results 
the following static equation: 

h(x,y) = -x x + x 2 + x 3 c x3 + 10" y - 

1Q y-PK w _ o 

2 + w PK 2-y 

C * 3 = 1 + 10 p *2-y + IQPK-L+PKi-ly 

PK X = - log 10 Ka-t 
PK 2 = - log 10 Ka 2 
The dynamic equations are given by: 

On 



dt V 
d± 2 _ <fi 



x 1 ) + f(w 2t 



y (W 12 -X 2 )+y (W 22 - X 2 ) + y (« 2 

dx 3 



dt 

Where 



t/2 



,. v 13 - *3J t — (w 23 - X 3 ) +y^;, 



x t ) 

X 2 ) 

Q3 (« 3 -X 3 ) 



(4) 
(5) 
(6) 



V: Volume of the mixing tank, ml 

Kw: Dissociation constant of water, 10~ 14 

Ka L : ith dissociation constant of acid 

Wj : Concentration of the ith species in the process 
stream, mol/1 

Wii: Concentration of the ith species in the acid stream, 
mol/1 

w 2i : Concentration of the ith species in the buffer 
stream, mol/1 

q t : Flow rate of acid, buffer and base stream in 
simulation, ml/s 

OC;: Concentration of the ith species in the titrating 
stream, mol/1 

Xi : Reaction invariant of ith species, mol/1 
y: Process variable, pH 

u: Flow rate of the titrating stream, ml/min or ml/s 
3. Simulation Results 

The closed loop control system was solved using 
Levenberg-Marquardt's optimization approach with sampling 
time of 0.001 s. The simulation method combines SIMULINK 
module for pH neutralization model and M-file for LMA 
approach. A list of M-file programs used in the paper is provided 
in Appendix 1 and 2. Figure (1 and 2) show the responses of the 
pH neutralization obtained with change in pH set point. The 
optimal gains of PID controller are calculated to minimize the 
error function which described in equation (2). Also the values 
of gains of PID controller are plotted. The response of pH value 
tends to set point value with minimum steady state error and the 
values of gains tend to minimum values once the error reaches to 
zero value. Figure (3) show the change in set point from pH = 6 
to 7, noticed that the spike value of manipulated variable (q 3 ) and 
the values of gains (Kp and Ki) tend to increasing while the gain 
(Kd) tends to decreasing. In other hand, figure (4) show the 
change in set point from 7 to 6, the value of manipulated variable 
(q 3 ) and the value of gains (Kp, Ki and Kd) tend to decreasing. 
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Fig. 1. Simulation results of PID controller and values of Fig. 2. Simulation results of PID controller and values of 
parameters parameters 
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Fig. 3. Simulation results of PID controller and values of 
parameters 
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Fig. 4. Simulation results of PID controller and values of 
parameters 
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CONCLUSION 

The paper presents an application of the Levenberg- 
Marquardt Algorithm (LMA) to optimization of parameters (Kp, 
Ki, Kd) of the PID controller structure according to minimum of 
integral square error. The simulated results were obtained of 
parameters by means of computer program implemented in 
Matlab software. As an example, the optimization of parameters 
of PID controllers with reference to a ph Neutralization process 
was presented. 
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Appendix (1) 

%Determination of PID controller prameters using "Levenberg- 
Marquardt %Algorithm " 
%Step 1: Write an M-file tracklsq.m. 
function [Kp,Ki,Kd] = phpidL2015k 
ph2015 % Load the simulink model 
k0 = [0 0 0]; % Set initial values of parametrs 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%%%%% 

%%%Create or edit optimization options structure 

options = optimset(' Algorithm', 'levenberg- 

marquardt', 'Display', 'iter',... 

'TolX',le-5, 'TolFun ',le-9, 'TolCon \le-6); 

k = lsqnonlin(@tracklsq, kO, [], [], options); 

Kp=k(l ); Ki = k(2 );Kd = k(3 ); 

function F = tracklsq(k) 

Kp=k(l); 

Ki = k(2); 

Kd = k(3); 
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%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%%%% 

% Choose solver and set model workspace to this function 
%Step 2: Invoke optimization routine, 
simopt = simset('solver', 'ode5', 'SrcWorkspace', 'Current'); 
[t,x,yl,y2,y3] = sim('ph2015',[0 400], simopt); 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%%%% 

%F = set point - actual value; % Compute error value 

F=y2-yl; 

end 

% 

% Put variables back in the base workspace 
Kp=k(l) 
Ki = k(2) 
Kd = k(3) 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%%% 

%plot( toutjout, 'r',tout,yset, '-b '); 
% Set axes and labels. 

%axis([0 30 -1.6 1.6]); xlabel('Time'); ylabel( 'Amplitude'); 
elf reset 

subplot(2,l,l ),plot( t,yl, 'r ', t,y2, '-.b', 'Line Width ',2); 
title( 'levenberg-marquardt') 

legend({'pH', 'pH setpoint'}, 'Location', 'SE', 'FontSize',8') 
xlabel( 'Time (min) ', 'FontSize ',12') 
ylabel( 'pH', 'FontSize ',12) 
grid on 

subplot(2,l,2),plot(t,y3, 'r', 'LineWidth' ,2); 
xlabel('Time (min)', 'FontSize', 12') 
ylabel( 'q3 (mils)', 'FontSize ', 12 ') 
grid on 
end 

Appendix (2) 

% Read Mat. Files of PID gains 

% ki(First number, second number)first number points to raw number 
% second number points to col. number 

% 

load ki.mat; %# assume this contains a matrix called ki 
load kp.mat 
load kd.mat 

for i=l:l:51 %# numbers of columns 

x(i) = ki(l,i); %# numbers from all rows, column 1, into X 

xl(i)=kp(l,i); 

x2(i)=kd(l,i); 

y(i)= ki(2,i); %# numbers from all rows, column 2, into Y 

yl(i)= kp(2,i); 

y2(i)=kd(2,i); 

subplot(3,l,l), (plot(x,y, 'o ') ); 
title( 'PID Parameters ', 'FontSize ',12') 
ylabel( 'Ki ', 'FontSize ',12) 
subplot(3,l,2),(plot(xl,yl, '+')); 
ylabel( 'Kp ', 'FontSize ',12) 
subplotf 3, 1,3 ), (plot(x2,y2, '* ') ); 
xlabel('Time (min)', 'FontSize', 12') 
ylabel('Kd', 'FontSize', 12) 

end 
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